#This to obtain:
#table_A11.tex
#table_A14.tex

library("broom")
library("dplyr")
library("grid")
library("tidyr")
library("sandwich")
library("lfe")
library("broom")
library("dplyr")
library("tidyr")
library("grid")
library("stargazer")
library("lmtest")
library("multiwayvcov")
library("ggplot2")
library("glue")
library("spdep")
library("spatialreg")
require("pscl")
require("boot")
require("MASS")
library("readxl")


#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")


########################################################################
#Function for returning results that include controls for SP. Autocorr.#
########################################################################


spatial_regression_glmnb <- function(data_sample, y, x, unique_observation_id) {
  clean_x = clean_list(x, "quad")
  clean_x = clean_list(clean_x, "bfe")
  
  
  specification = as.formula(paste(glue("{y}~"), 
                                   paste(clean_x, collapse = "+")))
  bw<-subset(data_sample, select =c(x, y, unique_observation_id))
  #Keeing track on row indices
  bw<-cbind(new_col_id = rownames(bw), bw)
  #Removing NAs
  bw<-bw[!sf::st_is_empty(bw), ]
  bw<-na.omit(bw) 
  #Step1: the function "poly2nb" builds a neighbours list based on
  #regions with contiguous boundaries, that is sharing one or more boundary point.
  list_nb=poly2nb(bw, row.names = dplyr::pull(bw, unique_observation_id))
  #Step2: Convert nb to listw type
  list_nb_w=tryCatch(nb2listw(list_nb, zero.policy=TRUE), error=function(err) NA)
  #Step3: Spatial Filtering
  reg1_res_sp<-tryCatch(SpatialFiltering(specification,
                                         data = bw, 
                                         nb = list_nb, 
                                         style="W", 
                                         ExactEV = F, 
                                         zero.policy = T,
                                         na.action=na.exclude),
                        error=function(err) NA)
  print("Done Spatial Filtering")
  #Step4: Include fitted results"
  bw_fitted<-tryCatch(cbind(bw, as.data.frame(fitted(reg1_res_sp)),
                            id = row.names(bw)), error=function(err) 
                            {print("NO NEIGHBORS"); data_sample})
  
  #Keeping track of indices
  #This will allow us to see exactly which observations were used in the regressions
  row.names(bw_fitted) <- bw_fitted$id
  bw_fitted<-subset(bw_fitted, select = -c(id))
  print(names(bw_fitted))
  list_fitted_ei<-names(bw_fitted[, grepl("vec", names(bw_fitted))])
  list_fitted_ei<-list_fitted_ei[list_fitted_ei != "Shape"]
  
  #Remove vec
  all_vars<-x
  specification_with_eigenv = as.formula(paste(glue("{y}~"), 
                                               paste(all_vars, collapse = "+")))
  
  
  print(specification_with_eigenv)
  model1 <- tryCatch(glm.nb(specification_with_eigenv, data = bw_fitted), error=function(err) NA)
  reg_moran<-moran.test(x=resid(model1),listw=list_nb_w, zero.policy = T)$statistic
  reg_moran_p<-moran.test(x=resid(model1),listw=list_nb_w, zero.policy = T)$p.value
  reg_moran_star<-paste(round(reg_moran, digits = 3), "",
                        ifelse(reg_moran_p<.01,"***",
                               ifelse(reg_moran_p<.05,"**",
                                      ifelse(reg_moran_p<.1,"*",
                                             ""))), "", "", sep = "")

  return(list(model1, reg_moran_star))
}



######################################
#Function for cleaning the list of FE#
#####################################

clean_list <- function(old_list, particle) {
  term <- paste(glue("^{particle}"))
  new_list = list() 
  for (i in old_list) {
    #print(i)
    temp_list = c()
    
    for (j in i){
      if (!grepl(term, j)){
        temp_list <- c(temp_list, j)
      }
    }
    l = length(new_list)
    new_list[[l+1]] <- temp_list
  }
  return(new_list)}


#Functions for coefficient plots
rerun_regression<-function(regression_name, string_dv, orig_data_for_clusters, cluster_name){
  #Process:
  #-takes already existing regression objects, associates them with an 
  #"official" label, calculates robust clustered SE and CI based on a cluster
  #that is in the original dataset
  #and puts the results in a tibble
  #-this will be used for ggplot later
  
  #Input:
  #-Regression model object after lm
  #-The label for the DV that you want displayed in ggplot
  #-original dataset that contains thhe names of the clusters
  
  #Output:
  #-A tibble with regression coefficients from the model that you ran with:
  #--robust SE
  #--robust CI
  
  #Recreate dataframe from the regression object
  remake_data<-regression_name$model
  #Take the string of the DV and IV
  dv<-regression_name$terms[[2]]
  ivs<-regression_name$terms[[3]]
  ivs_str<-paste(ivs)[2]
  dv_str<-paste(dv, "~", collapse = "")
  dv_ivs_str<-paste(dv_str, ivs_str, sep=" ")
  #Create specification and run model
  specification = as.formula(dv_ivs_str)
  model <- glm.nb(specification, data = remake_data)
  no_obs<-nobs(model)
  no_obs<-paste("Obs: ", as.character(no_obs), sep = "")
  #Calculate robust SE
  model_se<-se_robust_cluster(model, orig_data_for_clusters, cluster_name)
  #Create tibble
  model_tidy<-tidy(model)
  model_tidy$no_obs<-no_obs
  model_se_tidy<-tidy(model_se)
  names(model_se_tidy)[1] <- "term"
  names(model_se_tidy)[2] <- "robust.se"
  model_tidy2<-left_join(model_tidy, model_se_tidy, by = c("term"="term"))
  
  #Calculate CI
  df_ci<-as.data.frame(ci_robust_cluster(model, orig_data_for_clusters, cluster_name))
  df_ci <- cbind(term = rownames(df_ci), df_ci)
  names(df_ci)[2] <- "lb_95"
  names(df_ci)[3] <- "ub_95"
  
  #Merge CI to the dataframe
  model_tidy3<-left_join(model_tidy2, df_ci, by = c("term"="term"))
  
  #Only keep the treat variable
  model_tidy4<-subset(model_tidy3, term=="treat")
  model_tidy4$dv_official<-string_dv
  
  return(model_tidy4)
}


return_graphs_coef<-function(models_run, models_run_alias, orig_data_for_clusters, cluster_name){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" labels - models_run_alias, 
  #create cluster robust SE with name clusters taken from the original dataset - orig_data_for_clusters
  #keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  #-Original dataset name from which te function takes the cluster names
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    print(paste("Model", i, sep=" "))
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    model_tidy<-rerun_regression(x, name_x, orig_data_for_clusters, cluster_name)
    model_tidy$counter<-counter
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  irislabs1 <- rev(big_data$dv_official)
  irislabs2 <- rev(big_data$no_obs)
  mean_feature <- ggplot(data = big_data, aes(x=estimate, y=as.numeric(reorder(dv_official, desc(counter))))) +
    geom_point(size = 2) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0)+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs1),
                       labels = irislabs1,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    scale_x_continuous(name = "Estimate") +
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=14),
          axis.text.y = element_text(size=14),
          axis.title=element_text(size=16))  
  #return(big_data)
  return(mean_feature)
}

clean_col_labels<-function(column_labels){
  #Process:
  #-takes list of column labels for stargazer
  #and removes the unnecessary tags
  
  #Input:
  #-list of column labels for stargazer
  
  #Output:
  #-returns list of clean labels
  df<-list()
  length(df)<-length(column_labels)
  counter = 1
  for (i in column_labels){
    #print(i)
    j<-str_remove(i, '\\{')
    k<-str_remove(j, '\\}')
    l<-str_remove(k, '\\\\')
    m<-str_remove(l, 'shortstack')
    n<-gsub("\\\\", "\n", m)
    o<-gsub("\n\n", "\n", n)
    p<-gsub(",", "", o)
    df[counter] = p
    counter = counter +1
    df<-unlist(df)
  }
  return(df)
}

se_robust_cluster <- function(x, data_sample_shape, name_cluster){
  #Note: the model needs to include some variable whose values can be identified and
  #be put back in the original dataframe
  #Here this variable is dist_east_west_brd
  data_model<-x$model
  data_model$in_regression<-1
  data_model<-subset(data_model, select=c(in_regression))
  #Create index
  data_model <- tibble::rowid_to_column(data_model, "index")
  data2<-tibble::rowid_to_column(data_sample_shape, "index")
  #Merge by rowname or index
  data_merged <- left_join(data_model, data2, by=c("index" = "index"))
  #Remove any missing
  data_subset<-subset(data_merged, !is.na(in_regression))
  dat_clusters<-data_subset[[name_cluster]]
  res<-coeftest(x, vcov = cluster.vcov(x, dat_clusters, force_posdef=T))[, 2]
  return(res)}

ci_robust_cluster <- function(x, data_sample_shape, name_cluster){
  #Note: the model needs to include some variable whose values can be identified and
  #be put back in the original dataframe
  #Here this variable is dist_east_west_brd
  data_model<-x$model
  data_model$in_regression<-1
  data_model<-subset(data_model, select=c(in_regression))
  #Create index
  data_model <- tibble::rowid_to_column(data_model, "index")
  data2<-tibble::rowid_to_column(data_sample_shape, "index")
  #Merge by rowname or index
  data_merged <- left_join(data_model, data2, by=c("index" = "index"))
  #Remove any missing
  data_subset<-subset(data_merged, !is.na(in_regression))
  dat_clusters<-data_subset[[name_cluster]]
  res_ci<-coefci(x, vcov = cluster.vcov(x, dat_clusters, force_posdef=T))
  return(res_ci)}


#Reading Polygons
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

krajna <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3$ID_2

HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)


final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
final_sp$bosnia_border_mun
data<-subset(final_sp, bosnia_border_mun==0)


specifications <- list(c('treat', 
                         'POINT_X', 'POINT_Y', 
                         'zagreb_distance', 
                         'bfe1', 'bfe2', 'bfe3', 'bfe4', 'bfe5', 'bfe6', 'bfe7',
                         'quad1', 'quad2', 'quad3'))


#Model1
data$log_sum_best<-log(data$Sum_best + 1)
log_sum_best_list<-spatial_regression_glmnb(data, "log_sum_best", specifications[[1]], "ID_2")
m_log_sum_best<-log_sum_best_list[[1]]
log_sum_best_m<-log_sum_best_list[[2]]
summary(m_log_sum_best)

#Model2
sum_minor_dam_1_list<-spatial_regression_glmnb(data, "Sum_minor_dam_1", specifications[[1]], "ID_2")
m_sum_minor_dam_1<-sum_minor_dam_1_list[[1]]
sum_minor_dam_1_m<-sum_minor_dam_1_list[[2]]
summary(m_sum_minor_dam_1)

#Model3
sum_light_dam_constr_2_list<-spatial_regression_glmnb(data, "Sum_light_dam_constr_2", specifications[[1]], "ID_2")
m_sum_light_dam_constr_2<-sum_light_dam_constr_2_list[[1]]
sum_light_dam_constr_2_m<-sum_light_dam_constr_2_list[[2]]
summary(m_sum_light_dam_constr_2)

#Model4
sum_light_dam_port_3_list<-spatial_regression_glmnb(data, "Sum_light_dam_port_3", specifications[[1]], "ID_2")
m_sum_light_dam_port_3<-sum_light_dam_port_3_list[[1]]
sum_light_dam_port_3_m<-sum_light_dam_port_3_list[[2]]
summary(m_sum_light_dam_port_3)

#Model5
sum_serious_dam_4_list<-spatial_regression_glmnb(data, "Sum_serious_dam_4", specifications[[1]], "ID_2")
m_sum_serious_dam_4<-sum_serious_dam_4_list[[1]]
sum_serious_dam_4_m<-sum_serious_dam_4_list[[2]]
summary(m_sum_serious_dam_4)

#Model6
sum_build_part_destr_5_list<-spatial_regression_glmnb(data, "Sum_build_part_destr_5", specifications[[1]], "ID_2")
m_sum_build_part_destr_5<-sum_build_part_destr_5_list[[1]]
sum_build_part_destr_5_m<-sum_build_part_destr_5_list[[2]]
summary(m_sum_build_part_destr_5)

#Model7
sum_build_compl_dam_6_list<-spatial_regression_glmnb(data, "Sum_build_compl_dam_6", specifications[[1]], "ID_2")
m_sum_build_compl_dam_6<-sum_build_compl_dam_6_list[[1]]
sum_build_compl_dam_6_m<-sum_build_compl_dam_6_list[[2]]
summary(m_sum_build_compl_dam_6)

#Model8
sum_buildings_list<-spatial_regression_glmnb(data, "Sum_buildings", specifications[[1]], "ID_2")
m_sum_buildings<-sum_buildings_list[[1]]
sum_buildings_m<-sum_buildings_list[[2]]
summary(m_sum_buildings)

#Model9
sum_church_faith_places_list<-spatial_regression_glmnb(data, "Sum_church_faith_places", specifications[[1]], "ID_2")
m_sum_church_faith_places<-sum_church_faith_places_list[[1]]
sum_church_faith_places_m<-sum_church_faith_places_list[[2]]
summary(m_sum_church_faith_places)

#Model10
sum_schools_hotels_list<-spatial_regression_glmnb(data, "Sum_schools_hotels", specifications[[1]], "ID_2")
m_sum_schools_hotels<-sum_schools_hotels_list[[1]]
sum_schools_hotels_m<-sum_schools_hotels_list[[2]]
summary(m_sum_schools_hotels)

#Model11
sum_monuments_symbols_list<-spatial_regression_glmnb(data, "Sum_monuments_symbols", specifications[[1]], "ID_2")
m_sum_monuments_symbols<-sum_monuments_symbols_list[[1]]
sum_monuments_symbols_m<-sum_monuments_symbols_list[[2]]
summary(m_sum_monuments_symbols)


#######
#GRAPH#
#######

######################
#Step1: RELEVANT VARS#
######################

relevant_vars<-c("log_sum_best", 
                 "Sum_minor_dam_1", 
                 "Sum_light_dam_constr_2",
                 "Sum_light_dam_port_3", 
                 "Sum_serious_dam_4", 
                 "Sum_build_part_destr_5",
                 "Sum_build_compl_dam_6", 
                 "Sum_buildings", 
                 "Sum_church_faith_places",
                 "Sum_schools_hotels", 
                 "Sum_monuments_symbols")


#####################
#Step2: MEANS AND SD#
#####################

means<-list()
sds<-list()
for (i in relevant_vars){
  item_mean<-round(mean(data[[i]], na.rm=T), 3)
  item_sd<-round(sd(data[[i]], na.rm=T), 3)
  means[[length(means) + 1]] <- item_mean 
  sds[[length(sds) + 1]] <- item_sd 
}
means<-unlist(means)
sds<-unlist(sds)

####################
#Step3: LIST MODELS#
####################


list_models<-list(m_log_sum_best,
                  m_sum_minor_dam_1, 
                  m_sum_light_dam_constr_2,
                  m_sum_light_dam_port_3,
                  m_sum_serious_dam_4,
                  m_sum_build_part_destr_5,
                  m_sum_build_compl_dam_6, 
                  m_sum_buildings,
                  m_sum_church_faith_places,
                  m_sum_schools_hotels,
                  m_sum_monuments_symbols)

#####################
#Step3: COLUMN NAMES#
#####################


column_labels= c("\\shortstack{No. of\\\\ Deaths}", 
                 "\\shortstack{Light\\\\ Damage 1}",
                 "\\shortstack{Light\\\\ Damage 2}",
                 "\\shortstack{Light\\\\ Damage 3}", 
                 "\\shortstack{Serious\\\\ Damage 4}",
                 "\\shortstack{Partial\\\\ Destr. 5}",
                 "\\shortstack{Complete\\\\ Destr. 6}",
                 "\\shortstack{Buildings}", 
                 "\\shortstack{Churches}",
                 "\\shortstack{Schools}",
                 "\\shortstack{Monuments}")


yug_war<-stargazer(list_models, 
                   object.names = F,
                   keep = c("treat"),
                   se = list(se_robust_cluster(m_log_sum_best, data, "NAME_2"),
                             se_robust_cluster(m_sum_minor_dam_1, data, "NAME_2"),
                             se_robust_cluster(m_sum_light_dam_constr_2, data, "NAME_2"),
                             se_robust_cluster(m_sum_light_dam_port_3, data, "NAME_2"),
                             se_robust_cluster(m_sum_serious_dam_4, data, "NAME_2"),
                             se_robust_cluster(m_sum_build_part_destr_5, data, "NAME_2"),
                             se_robust_cluster(m_sum_build_compl_dam_6, data, "NAME_2"),
                             se_robust_cluster(m_sum_buildings, data, "NAME_2"),
                             se_robust_cluster(m_sum_church_faith_places, data, "NAME_2"),
                             se_robust_cluster(m_sum_schools_hotels, data, "NAME_2"),
                             se_robust_cluster(m_sum_monuments_symbols, data, "NAME_2")),
                   column.labels= column_labels,
                   covariate.labels=c("Military Territory"),
                   dep.var.labels.include = F,
                   omit.stat=c("LL","ser","f", "rsq", "theta"), no.space=TRUE, float=FALSE,
                   add.lines=list(c("Mean", means),
                                  c("SD", sds),
                                  c("Boundary FE", rep("Yes", 14)),
                                  c("Region FE", rep("Yes", 14))))


note_latex <- "\\multicolumn{12}{l} {\\parbox[t]{28cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and standard errors in parantheses from a negative binomial. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis is municipality (opcina). See codebook for data sources and how the variables were calculated.}}"
yug_war [grepl("Note", yug_war)] <- note_latex
cat(yug_war, file="./Dropbox/Legacies_Project/Paper/tables/table_A14.tex", sep="\n")



#####
#WW2#
#####

#Model11
log_ww2_bombs_no_list<-spatial_regression_glmnb(data, "log_ww2_bombs_no", specifications[[1]], "ID_2")
m_log_ww2_bombs_no<-log_ww2_bombs_no_list[[1]]
log_ww2_bombs_no_m<-log_ww2_bombs_no_list[[2]]
summary(m_log_ww2_bombs_no)

log_tons_TNT_list<-spatial_regression_glmnb(data, "log_tons_TNT", specifications[[1]], "ID_2")
m_log_tons_TNT<-log_tons_TNT_list[[1]]
log_tons_TNT_m<-log_tons_TNT_list[[2]]
summary(m_log_tons_TNT)

log_con_camp_no_list<-spatial_regression_glmnb(data, "log_con_camp_no", specifications[[1]], "ID_2")
m_log_con_camp_no<-log_con_camp_no_list[[1]]
log_con_camp_no_m<-log_con_camp_no_list[[2]]
summary(m_log_con_camp_no)

#log_con_camp_prisoners_list<-spatial_regression_glmnb(data, "log_con_camp_prisoners", specifications[[1]], "ID_2")
#m_log_con_camp_prisoners<-log_con_camp_prisoners_list[[1]]
#log_con_camp_prisoners_m<-log_con_camp_prisoners_list[[2]]
#summary(m_log_con_camp_prisoners)

m_log_con_camp_prisoners<-lm(log_con_camp_prisoners~treat + POINT_X + POINT_Y + zagreb_distance +
                              bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                              quad1 + quad2 + quad2, data = data)
summary(m_log_con_camp_prisoners)       


#log_con_camp_deaths_list<-spatial_regression_glmnb(data, "log_con_camp_deaths", specifications[[1]], "ID_2")
#m_log_con_camp_deaths<-log_con_camp_deaths_list[[1]]
#log_con_camp_deaths_m<-log_con_camp_deaths_list[[2]]
#summary(m_log_con_camp_deaths)

m_log_con_camp_deaths<-lm(log_con_camp_deaths~treat + POINT_X + POINT_Y + zagreb_distance +
                               bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                               quad1 + quad2 + quad2, data = data)
summary(m_log_con_camp_deaths)       



# log_massacre_death_no_list<-spatial_regression_glmnb(data, "log_massacre_death_no", specifications[[1]], "ID_2")
# m_log_massacre_death_no<-log_massacre_death_no_list[[1]]
# log_massacre_death_no_m<-log_massacre_death_no_list[[2]]
# summary(m_log_massacre_death_no)

m_log_massacre_death_no<-lm(log_massacre_death_no~treat + POINT_X + POINT_Y + zagreb_distance +
         bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
         quad1 + quad2 + quad2, data = data)
summary(m_log_massacre_death_no)       

m_log_massacre_death_no<- lm(log_massacre_death_no ~treat+
                                 POINT_X+POINT_Y+zagreb_distance+
                                 bfe1+bfe2+bfe3+bfe4+
                                 bfe5+bfe6+bfe7+
                               quad1+ quad2 + quad3,
                               data=data)

########
#GRAPH#
#######


######################
#Step1: RELEVANT VARS#
######################

relevant_vars<-c("log_ww2_bombs_no", 
                 "log_tons_TNT", 
                 "log_con_camp_no",
                 "log_con_camp_prisoners", 
                 "log_massacre_death_no")


#####################
#Step2: MEANS AND SD#
#####################


means<-list()
sds<-list()
for (i in relevant_vars){
  item_mean<-round(mean(data[[i]], na.rm=T), 3)
  item_sd<-round(sd(data[[i]], na.rm=T), 3)
  means[[length(means) + 1]] <- item_mean 
  sds[[length(sds) + 1]] <- item_sd 
}
means<-unlist(means)
sds<-unlist(sds)

####################
#Step3: LIST MODELS#
####################


list_models<-list(m_log_ww2_bombs_no, 
                  m_log_tons_TNT, 
                  m_log_con_camp_no,
                  m_log_con_camp_prisoners, 
                  m_log_massacre_death_no)


column_labels= c("\\shortstack{Log. No.\\\\ Bombs}",
                 "\\shortstack{Log. Tons\\\\ TNT}",
                 "\\shortstack{Log Concentration\\\\ Camp Count}",
                 "\\shortstack{Log Concentration\\\\ Camp prisoners}",
                 "\\shortstack{Log Massacre\\\\ Death Count}")



ww2 =stargazer(list_models, 
               object.names = F,
               model.names = F,
               keep = c("treat"),
               se = list(se_robust_cluster(m_log_ww2_bombs_no, data, "NAME_2"),
                         se_robust_cluster(m_log_tons_TNT, data, "NAME_2"),
                         se_robust_cluster(m_log_con_camp_no, data, "NAME_2"),
                         se_robust_cluster(m_log_con_camp_prisoners, data, "NAME_2"),
                         se_robust_cluster(m_log_massacre_death_no, data, "NAME_2")),
               column.labels= column_labels,
               covariate.labels=c("Military Territory"),
               dep.var.labels.include = F,
               omit.stat=c("LL","ser","f", "rsq", "theta"), no.space=TRUE, float=FALSE,
               add.lines=list(c("Mean", means),
                              c("SD", sds),
                              c("Boundary FE", rep("Yes", 14)),
                              c("Region FE", rep("Yes", 14))))
note_latex <- "\\multicolumn{6}{l} {\\parbox[t]{18cm}{ \\footnotesize \\textit{Notes:} 
Results are coefficients and standard errors from negative binomial models.
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis is municipality (opcina). See codebook for data sources and how the variables were calculated.}}"
ww2 [grepl("Note", ww2)] <- note_latex
cat(ww2, file="./Dropbox/Legacies_Project/Paper/tables/table_A11.tex", sep="\n")

